%% PRACTICAL: MANUAL PREPROCESSING 
% 
% In this practical we will take an approach which is run step-by-step and 
% requires manual intervention. This will go through the following steps:
% 
%     1) Conversion of data into SPM format
%     2) Downsampling
%     3) High pass filtering
%     4) Removing bad time periods using OSLview
%     5) Manual AFRICA denoising
%     6) Epoching
%     7) Semi-automated outlier trial rejection (using a Fieldtrip tool) 
%
% We will work with a single subject's data from an emotional faces task 
% (data courtesy of Susie Murphy). These can be downloaded from:
% 
% <http://www.fmrib.ox.ac.uk/~woolrich/faces_subject1_data.tar.gz>
%
% Note that this contains the fif file:
% 
% fifs/sub1_face_sss.fif
%
% Which has already been SSS Maxfiltered and downsampled to 250 Hz, and 
% which we will use as the input into this analysis.
%
% MWW & APB 2014


%% SETUP THE MATLAB PATHS
%
% Sets the Matlab paths to include OSL. Change these paths so that they 
% correspond to the setup on your computer. You will also need to ensure 
% that fieldtrip and spm are not in your matlab path (as they are included
% within OSL).
    
osldir = '/Users/abaker/Dropbox/Code/osl_latest/osl_svn/osl';    

addpath(osldir);
osl_startup(osldir);


%% SPECIFY DIRECTORIES FOR THIS ANALYSIS
%
% Change the datadir and workingdir variables to correspond to the correct
% directories.

% This is the directory containing the downloaded experimental data
datadir = '/Users/abaker/Scratch/faces_subject1_data';
    
% This is the directory in which the analysis files will be stored
workingdir = datadir;

if ~isdir(workingdir)
    mkdir(workingdir);
end


%% SET UP THE LIST OF SUBJECTS AND THEIR STRUCTURAL SCANS FOR THE ANALYSIS
%
% Specify a list of the fif files, structual files (not applicable for this
% practical) and SPM files (which will be created). It is important to make
% sure that the order of these lists is consistent across sessions. 
% Note that here we only have 1 subject, but more generally there would be
% more than one, e.g.:

% fif_files{1}=[testdir '/fifs/sub1_face_sss.fif']; 
% fif_files{2}=[testdir '/fifs/sub2_face_sss.fif']; 
% etc...

% spm_files{1} = [workingdir '/sub1_face_sss.mat'];
% spm_files{2} = [workingdir '/sub2_face_sss.mat'];
% etc...

clear fif_files spm_files_basenames;

% List of the existing fif files for each subject
fif_files{1} = [datadir '/fifs/sub1_face_sss.fif']; 

% List of SPM MEEG object files to be created
spm_files_basenames{1} = 'sub1_face_sss.mat';

%% CONVERT FROM FIF TO AN SPM MEEG OBJECT:
%
% The fif file that we are working with is sub1_face_sss.fif. This has
% already been maxfiltered for you and downsampled to 250Hz.
% 
% This converts the data into an SPM MEEG object, and displays a histogram 
% plot showing the number of events detected for each code on the trigger 
% channel. The codes used on the trigger channel for this experiment were: 
%
%  1 = Neutral face
%  2 = Happy face
%  3 = Fearful face
%  4 = Motorbike
% 18 = Introduction screen
% 11 = Break between blocks
% 19 = Midway break
% 12 = Green fixation cross (response trials)
% 13 = Red fixation cross (following green on response trials)
% 14 = Red fixation cross (non-response trials)
%
% Note that there should be 120 motorbike trials, and 80 of each of the 
% face conditions. Check that the histogram plot corresponds to these trial 
% numbers.

spm_files = cell(size(fif_files));
for subnum = 1:length(fif_files), % iterates over subjects
    spm_files{subnum}=[workingdir '/' spm_files_basenames{subnum}];
end

if(~isempty(fif_files))
    
    S = [];

    for subnum = 1:length(fif_files), % loops over subjects
        
        S.fif_file = fif_files{subnum};
        S.spm_file = spm_files{subnum};       
        S.trigger_channel_mask = '0000000000111111'; % binary mask to use on the trigger channel
        
        % The conversion to SPM will show a histogram of the event codes
        % and correspond to those listed below in the epoching section
        [D, spm_files{subnum}] = osl_convert_script(S);
    end
end

%%%
% Note that this spmfile is the output from the conversion:
spm_files{1}


%% LOAD THE SPM MEEG OBJECT
%
% This will display the summary information about the SPM MEEG object. 
% Note that at this stage in the pipeline the data is continuous (with only 
% 1 "trial"), with 232,000 time points sampled at 250 Hz. We will epoch the 
% data later.

% Set new SPM MEEG object filenames to be used in following steps
for subnum = 1:length(spm_files), % iterates over subjects
    spm_files{subnum}=[workingdir '/' spm_files_basenames{subnum}];
end

% load in the SPM MEEG object
subnum = 1;
D = spm_eeg_load(spm_files{subnum});

% View the SPM MEEG object
display(D);

%% Remove Jumps (discontinuities)
%
% Any discontinuities in the data will cause large artefacts when passed
% through digital filters. They may also have associated ringing artefacts
% from the hardware filters in the scanner. 
%
% You can look at your data with 
% oslview(D);
% 
% Artefacts can be removed with 
% osl_remove_jumps(D);

%% DOWNSAMPLE
%
% This data has already been downsampled to 250 Hz when the Maxfilter was 
% run on it. But we will now do further downsampling as this helps to speed 
% things up even more, and we do not need information at high frequency in 
% this particular analysis. [Note that doing downsampling here is 
% particularly necessary if movement compensation has been used when 
% running Maxfilter, as this stops you from doing downsampling as part of 
% the Maxfilter call.]

% Set new SPM MEEG object filenames to be used in following steps
for subnum = 1:length(spm_files), % iterates over subjects
    spm_files{subnum}=[workingdir '/' spm_files_basenames{subnum}];
end

% Apply downsampling to each subject. Note that this will create new SPM
% files with the prefix 'd' (e.g. 'dsub1_face_sss.mat')
for subnum = 1:length(spm_files)
    S = [];
    S.D = spm_files{subnum};
    S.fsample_new = 150; % in Hz
    D = spm_eeg_downsample(S);    
end
close all

% Set new SPM MEEG object filenames to be used in following steps
for subnum = 1:length(spm_files), % iterates over subjects
    spm_files{subnum}=[workingdir '/d' spm_files_basenames{subnum}];
end

% Load in the SPM MEEG object
subnum = 1;
D = spm_eeg_load(spm_files{subnum});

% Look at the SPM object. Note that it is continuous data, with 139200 time
% points at 150 Hz. We will epoch the data later.
D


%% HIGH PASS FILTERING
%
% High-pass filter the data above 0.1 Hz to remove baseline drifts:

% Set new SPM MEEG object filenames to be used in following steps
for subnum = 1:length(spm_files), % iterates over subjects
    spm_files{subnum}=[workingdir '/d' spm_files_basenames{subnum}];
end

% Apply filtering to each subject. Note that this will create new SPM
% files with the prefix 'd' (e.g. 'fdsub1_face_sss.mat')
S              = [];
S.D            = spm_files{1};
S.filter.band  = 'high';
S.filter.PHz   = 0.1;

D=spm_eeg_filter_v2(S);


%% OSLVIEW
%
% Note that there are some large artefacts. Use the OSLview functionality
% to remove the bad epochs (see
% <https://sites.google.com/site/ohbaosl/preprocessing/oslview>)

% Set new SPM MEEG object filenames to be used in following steps
for subnum = 1:length(spm_files), % iterates over subjects
    spm_files{subnum}=[workingdir '/fd' spm_files_basenames{subnum}];
end

%%%
% *Load OSLview:*
% This data has some bad artefacts in. Mark the epochs at around 325s,
% 380s and 600s as bad. Plus mark the really bad artefacts at the end of
% the experiment from about 650 secs to the end. This will mean that we are 
% not using about half of the data. But with such bad artefacts this is the
% best we can do. We can still obtain good results with what remains.

for subnum = 1:length(spm_files)
    D = spm_eeg_load(spm_files{subnum});
    D = oslview(D);
    waitfor(gcf) % wait until current OSLview window has been closed
end


%% AFRICA WITH MANUAL COMPONENT REJECTION
%
% Run AFRICA ICA denoising. AFRICA uses independent component analysis to
% decompose sensor data into a set of maximally independent time courses.
% Using this framework, sources of interference such as eye-blinks, ECG
% artefacts and mains noise can be identified and removed from the data. 
%
% In this practical we will use manual artefact rejection by looking at the
% time courses and sensor topographies of each component and rejecting
% those that correlate with EOG and ECG measurements.
%
% The user interface displays the time course, power spectrum and sensor 
% topography for each component. These components are sorted based on one
% of a number of metrics, which you can toggle using the dropdown menu.


% Set new SPM MEEG object filenames to be used in following steps
for subnum = 1:length(spm_files), % iterates over subjects
    spm_files{subnum}=[workingdir '/fd' spm_files_basenames{subnum}];
end

%%%
% *AFRICA with manual artefact rejection:* Scroll through components using 
% the cursor keys. Identify the two components that correlate with the EOG 
% and ECG measurements and mark them for rejection using the red cross.
% Close the window when finished to save your results.
for subnum = 1:length(spm_files)
    
    [dirname,filename] = fileparts(spm_files{subnum});
    
    S = [];
    S.fname           = fullfile(dirname, filename);
    S.logfile         = 1;
    S.ica_file        = fullfile(dirname,[filename '_africa']);
    S.used_maxfilter  = 1;
    S.ident.func      = @identify_artefactual_components_manual;
    S.to_do           = [1 1 1]; 
    S.ident.artefact_chans = {'EOG','ECG'};
    
    osl_africa(S);
    
end

%% EPOCHING
%
% Does preliminary epoching for the purpose of finding outliers
% This is not the final epoching. Instead this sets up the epoch
% definitions, and performs a temporary epoching for the purpose of doing
% semi-automated outlier trial rejection (before running the fully 
% automated OAT).
%
% The epoch definitions and the continuous data will be kept and
% passed into OAT. This is so that things like temporal filtering (which is
% dones as part of OAT) can be done on the continuous data, before the data
% is epoched inside OAT.
%
% Note that this will also remove those trials that overlap with the bad
% epochs identified using OSLview. 
%
% Here the epochs are set to be from -1000ms to +2000ms relative to the 
% triggers in the MEG data. We also specify the trigger values for each of 
% the 4 epoch types of interest (motorcycle images, neutral faces, fearful 
% faces, happy faces). 
%
% The codes used on the trigger channel for this experiment were:
%  1 = Neutral face
%  2 = Happy face
%  3 = Fearful face
%  4 = Motorbike
% 18 = Introduction screen
% 11 = Break between blocks
% 19 = Midway break
% 12 = Green fixation cross (response trials)
% 13 = Red fixation cross (following green on response trials)
% 14 = Red fixation cross (non-response trials)
%
% Note that we're only interested in the first 4 event codes listed here for today's workshop. 

% Set new SPM MEEG object filenames to be used in following steps
for subnum = 1:length(spm_files), % iterates over subjects
    spm_files{subnum}=[workingdir '/Afd' spm_files_basenames{subnum}];
end

% Apply epoching to each subject. Note that this will create new SPM
% files with the prefix 'd' (e.g. 'efdsub1_face_sss.mat')
for subnum=1:length(spm_files)

    % Load dataset
    D_continuous = spm_eeg_load(spm_files{subnum});
    
    % Define the trials we want from the event information
    S = [];
    S.pretrig  = -1000; % epoch start in ms
    S.posttrig =  2000; % epoch end in ms  
    
    S.trialdef(1).conditionlabel = 'Neutral face';
    S.trialdef(1).eventtype      = 'STI101_down';
    S.trialdef(1).eventvalue     = 1;
    S.trialdef(2).conditionlabel = 'Happy face';
    S.trialdef(2).eventtype      = 'STI101_down';
    S.trialdef(2).eventvalue     = 2;
    S.trialdef(3).conditionlabel = 'Fearful face';
    S.trialdef(3).eventtype      = 'STI101_down';
    S.trialdef(3).eventvalue     = 3;
    S.trialdef(4).conditionlabel = 'Motorbike';
    S.trialdef(4).eventtype      = 'STI101_down';
    S.trialdef(4).eventvalue     = 4;
    
    S.reviewtrials      = 0;
    S.save              = 0;
    S.epochinfo.padding = 0;
    S.event             = D_continuous.events;
    S.fsample           = D_continuous.fsample;
    S.timeonset         = D_continuous.timeonset;
    
    [epochinfo.trl, epochinfo.conditionlabels] = spm_eeg_definetrial(S);        
    
    
    % Do epoching
    S2 = [];
    S2.D = D_continuous;     
    S2.epochinfo=epochinfo;
    [D,good_trial_start_times] = osl_epoch(S2);
            
end


%% EXAMINING THE EPOCHED DATA
%
% Now we'll take a look at our preprocessed data so far. Note that this is 
% now EPOCHED data, with:
% 4 conditions
% 323 channels
% 451 samples per trial
% 360 trials

% Set new SPM MEEG object filenames to be used in following steps
for subnum = 1:length(spm_files), % iterates over subjects
    spm_files{subnum}=[workingdir '/eAfd' spm_files_basenames{subnum}];
end

% Load in the SPM MEEG object
subnum = 1;             
D = spm_eeg_load(spm_files{subnum});

% Look at the SPM object. Note that this is now EPOCHED data.
D

% Display a list of trial types:
D.condlist

% Display time points (in seconds) per trial
D.time

% Identify trials of a certain type using the pickcondition function. E.g.:
motorbike_trls = D.pickconditions('Motorbike')

%%%
% Identify channels of certain types using the meegchannels function. E.g. 
% identify the channel indices for the planar gradiometers (Note that you 
% can use 'MEGMAG' to get the gradiometers, and D.chantype gives you a list 
% of all channel types by index).
planars = D.meegchannels('MEGPLANAR')

%%%
% We can access the actual MEG data using the syntax:
% D(channels, samples, trials). E.g. plot a figure showing all the trials
% for the motorbike condition in the 135th MEGPLANAR channel. Note that the 
% squeeze function is needed to remove single dimensions for passing to the 
% plot function, and D.time is used to return the time labels of the within 
% trial time points in seconds. 
figure; plot(D.time,squeeze(D(planars(135),:,motorbike_trls))); xlabel('time (seconds)');

%%%
% We can average over all the motorbike trials to get a rudimentary ERF
% (Event-Related Field):
figure; plot(D.time,squeeze(mean(D(planars(135),:,motorbike_trls),3))); xlabel('time (seconds)');

% Although we should bear in mind that this data is averaging over all data
% including noisy data segments, channels and trials. To do better than 
% this we need to perform outlier rejection.



%% VISUAL ARTEFACT REJECTION 
% This runs a Fieldtrip interactive tool.
%
% - Pass over the first interactive figure as it is the EOG channel - so just
% press the "quit" button.
%
% This will bring up another interactive figure which will show the
% magnetometers. You can choose the metric to display - it is best to stick
% to the default, which is variance. This metric is then displayed for 
% the different trials (bottom left), the different channels (top
% right), and for the combination of the two (top left). You need to use
% this information to identify those trials and channels with high variance
% and remove them.  
%
% 1) Remove the worst channel (with highest variance) by drawing a box
%    around it in the top right plot with the mouse. 
% 
% 2) Now remove the trials with high variance by drawing a box
%    around them in the bottom left plot.
%
% 3) Repeat this until you are happy that there are no more outliers.
%
% 4) Press "quit" and repeat the process for the gradiometers.

% Set new SPM MEEG object filenames to be used in following steps
for subnum = 1:length(spm_files), % iterates over subjects
    spm_files{subnum}=[workingdir '/eAfd' spm_files_basenames{subnum}];
end

%%%
% Run the visual artefact rejection tool:
for subnum=1:length(spm_files),
    S = [];
    S.D = spm_files{subnum};
    S.time_range = [-0.2 0.4];
    D = osl_rejectvisual(S);
end


%% EXAMINE THE CLEANED EPOCHED DATA
%
% We can now repeat the average over all the motorbike trials with the bad 
% trials removed to get a rudimentary ERF (Event-Related Field). 

% Set new SPM MEEG object filenames to be used in following steps
for subnum = 1:length(spm_files), % iterates over subjects
    spm_files{subnum}=[workingdir '/eAfd' spm_files_basenames{subnum}];
end

% load in SPM MEEG object
subnum = 1;             
D = spm_eeg_load(spm_files{subnum});

% List the marked bad channels
D.badchannels

% List the marked bad trials
D.reject

%%%
% Identify the motorbike trials. Note that pickconditions only ever 
% includes good trials.
motorbike_trls = D.pickconditions('Motorbike')

%%%
% Plot a cleaned rudimentary ERF
figure; plot(D.time,squeeze(mean(D(planars(135),:,motorbike_trls),3))); xlabel('time (seconds)');

%%%
% In practise you can use tools like OAT to do these ERF and cleverer 
% analyses over all sensors, and in source space.




